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Abstract 


In practical cases when a plate is subjected to dynamic loading, the geometric and 
material linearity may not be preserved. This study is aimed at investigating the effect of 
geometric non-linearity in the dynamic response. As there are no closed form solutions 
available owing to the complexity of the problem, numerical methods are used to solve the 
problem. 

For this computational woiic an efficient Finite Element code is developed using 
the “robust triangular elementf’ developed by Zien]dewiz[1991]. “Mindhn plate theory” in 
the form of mixed formulation is used for modeling plate bending. Inertia effect of the 
element is represented by consistent mass matrix. When the displacements are finite the 
effect of geometric non-linearity is not small and it has to be taken into consideration. 

Due to geometric non-linearity the structure stiffens and it increases with displacement 
Since the dynamic loading involves finite displacements, the effect of geometric non- 
linearity is analyzed in the present work. 

The computer code is first used to solve a standard linear static problem and after 
its validation it was used to solve a standard problem considering geometric non-linearity, 
and then extended to solve the problem involving crack. 



CHAPTER 1 


INTRODUCTION 


1.1 INTRODUCTION 


In the analysis of many structural problems it is implicitly assumed that both 
displacements and strains developed are small. In practical terms it means that 
geometry of the elements remain basically unchanged during the loading process and 
first order infinitesimal linear strain approximations can be used. It also has the 
advantage of the behaviours of a plate under in-plane load and flexural load are 
independent of each other. 

In practice such assumptions fail frequently even if strains may be small and 
elastic limits of ordinary structural materials are not exceeded. If accurate analysis is 
to be done and failure models are to be developed geometric non-linearity have to be 
considered in structures. 

For instance stress due to membrane action usually neglected in plate flexure 
may cause considerable decrease of displacements as compared with linear solution 
even though displacements are still quite small. Further, when an in-plane 
compressive load reaches certain value, deflections increase more rapidly than 
predicted by linear solution and indeed a state may be attained where load carrying 
capacity decreases with continuing deformation. 
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The applications of such an analysis clearly is of considerable importance in 
aerospace applications, radio telescopes, cooling towers, and other relatively slender 
structures, and crack propagation in structures. 

Many structures are subjected to dynamic loading. These loads may be 
created by moving vehicles, wind gusts, seismic disturbances, unbalanced machines, 
wave impacts, blasts. While solving these problems often the time variation is 
neglected and they are treated as static or quasi-static problems. However, if the time 
taken by the load to reach the maximum value is significantly less than fundamental 
period of the structure, then the inertia effect play an important role and the problem 
is termed as dynamic. 

Geometric non-linearity effect is of significant value when displacements are 
large or strains are finite or both. Since many of dynamic loads are large giving rise 
to large displacements, geometric non-linearity has to betaken into consideration while 
solving such problems. The presence of cracks in plate structures subjected to 
suddenly applied loads like explosions makes study of dynamic fracture with geometric 
non-linearity necessary for developing sound design methodologies aimed at assuring 
the integrity of the structure and safety of the personnel. For many of the dynamic 
problems it is impossible to find out closed form solutions especially if they involve 
fracture. Thus, it often becomes mandatory to analyze dynamic crack propagation 
problems with geometric non-linearity, in finite solids using numerical methods. 
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1.2 LITERATURE SURVEY 


The thin plate theory developed by Kirchhoff is modified for thick plates by 
Reissner and Mindilin by relaxing the assumption that the normals to the mid plane 
remain normal. This was done since the shear force terms cannot be neglected in 
many loading support conditions. Reissner and Mindlin theories can be extended to 
thick plates and thus known as Reissner-Mindilin formulation. Regarding Finite 
Element Method difference between thin plate and thick plate analysis is that in thick 
plate formulation instead of one approximate function for w (for calculation of rotational 
angle) three independent functions for deflection w and rotations 8^, e„ are taken. In 
thin plates it is possible to represent the state of deformation by one quantity deflection 
w and thus its interpolation function is of a minimum continuity is needed, while in 
the thick plates the interpolation functions for w along with rotation and shear degree 
of freedom need to be of C° continuity. 

While analytical solutions remain difficult for thick plates, the advent of Finite 
Element Method made thick plate theory simpler to implement. It is more convenient 
to start with thick plate formulation and can be used to thin plates if necessary, by 
imposition of certain relation between degrees of freedom [Zienkiewicz and 
Taylor, 1991]. During formulation, if shear forces are eliminated, the formulation is 
known as irreducible formulation and if shear forces are retained it is known as hybrid 
formulation. The latter formulation gives more degree of freedom and so it is preferred 
in most of the cases [Zienkiewicz and Taylor, 1991]. 
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other important factor is the selection of suitable element to perform 
satisfactorily for both thick and thin plates. While many elements are attempted to 
avoid problems of locking and singularity, the results by the use of ’robust triangular 
element’ are found good [Zienkiewicz and Lefebvre,1988]. 

Generally in many structural problems it is assumed that both displacements 
and strains developed are small. This has the advantage of the behaviour of plate 
under in-plane and flexural loads independent of each other. The other important 
advantage is that the Finite Element analysis is far simpler and computationally faster 
as original configuration remains unchanged. In practice such assumptions fail 
frequently when strains are finite or deformations are large or both and elastic limits 
are not reached. For better designing ’Geometric non-linearity’ should be considered. 
The application of geometric non-linearity has considerable importance in aerospace 
engineering, design of radio telescopes, cooling towers and many other[Zienkiewic 2 
and Taylor ,1991]. 

I 

In recent years a need has developed in industry to obtain non-linear analysis 
of structures subjected to static and dynamic loads. The earliest non-linear finite 
element analysis were essentially based on extensions of linear analysis and have 
been developed for specific applications [Oden, 1972]. 

Closed form analytical solutions are available for only a limited range of 
structures and often involve gross simplifications. For most practical problems non- 
linear analysis can be carried out only by numerical techniques [Mondkar and 


4 



Poweli.1977]. 


Basically two different approaches have been pursued in incremental non-linear 
finite element analysis. In the first static and kinematic variables are referred to 
updated configuration in each load step called as Eulerian or moving coordinate 
formulation. In second approach which is generally called Lagrangian formulation ail 
static and kinematic variables are referred to initial configuration [Bathe, 1975]. 

On comparison of Eulerian and Lagrangian formulation for a simple truss 
structure a maximum difference of about 25 per cent in displacements was found 
[Yamada,1972]. It was found that Lagrangian formulation is more general and 
computational more efficient. 

Of the various formulations of non-linearity some are general and some are 
restricted to account for material non-linearity only or for large displacements. Limited 
results have been obtained in dynamic non-'linear analysis involving large 
displacements and finite strains. In dynamic analysis numerical time integration of finite 
element equations of motions is required .Extensive research is devoted towards the 
development of stable and accurate integration scheme. [Bathe, Ramm, Wilson, 1975]. 

One of the important areas of concern is the iterative scheme to be used in 
non-linear problems. Physical insight into nature of the problem and usually small step 
incremental approaches are essential to obtain physically significant answers. In 
general, non-linear response will be computed by a combination of step by step and 
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iterative procedures. Depending on the degree of non-linearity equilibrium iterations 
may are may not be needed, and tangent stiffness may or may not be formulated in 
every step. Because of large computational effort typically required for decomposition 
of tangent stiffness matrix, it is desirable to seek a strategy in which the number of 
stiffness matrix reformulations are minimized [Modkar and Powell, 1977]. 

Of the many iterative schemes to solve the non-linear problems Newton- 
Raphson method is found favourable in terms of ease in mathematical formulation, 
indeed it is the only process in which convergence is quadratic. [Simo and 
Taylor, 1986]. The Newton-Raphson process despite its rapid convergence can be 
expensive and inconvenient since stiffness matrix has to be formed and factorized for 
each integration step, some times symmetric matrix becomes non symmetric. These 
drawbacks can be overcome by modified Newton-Raphson method, but this method 
suffers from a draw back of converging slowly [Zienkiewicz and Taylor, 1991]. The 
convergence rate of modified Newton-Raphson method can be improved by selective 
relaxation iterative solution technique but this is highly problem-dependent. It will be 
very difficult to use this modification when both tension and compression are present 
[Carter Wellford and Halissen,1981]. 

It was only after the year 1960 that the stress problem of cracked plates are 
examined from a higher order plate bending theory taking shear deformation into 
account. The ealiest attempt [Knowles, 1960] was concerned with the problem of 
infanite plate under uniform uniaxial bending far from the crack and was restricted to 
the case of vanishingly small plate thickness. 
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The predictions of classical theory [Williams, 1975] are known to be erroneous 
with respect to the angular distribution of stress around the crack tip. Results reported 
[Hartranft,1968] show that the stress intensity factor varies very rapidly in the h/a 
range of 0.0 to 0.25. In paper [Murty, 1981] the fmiteness of the plate was taken into 
account and the analysis here was based on differential approach. Shear deformation 
was taken into account through the use of Reissner theory. No restriction was put on 
the h/a ratio so that thickness effect could automatically be taken into account. 

Numerical solutions based on the finite element method have been obtained for 
a number of elastodynamic crack growth problems in two dimensions. Nishioka and 
Atiuri [1986] gave a very elaborate information on analysis of dynamic fracture using 
finite element method. The most common ways to deal with the crack-tip region are 
to simulate crack growth through gradual release of element nodal forces or 
imbedding a moving element in which the interpolation functions are determined by 
the continuum near tip fields at the crack-tip in the mesh, and energy path integral 
considerations. 

Research at the Defence Research Establishment, Suffield on vulnerability of 
naval structures has been directed to study the dynamic response of unstiffened and 
stiffened panels subjected to air blast shock waves. Numerical modelling using the 
finite element code ADINA has been contributed significantly to this research and has 
provided benchmark solutions. It involves a thin four noded shell element [Houlston 
and Slater, 1991] and use of isotropic elasto-plastic bilinear hardening material. The 
study highlights the importance of nonlinear analysis and the detailed approximations 
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to the actual boundary conditions. 

A detailed analysis of the various time integration schemes from the point of 
view of accuracy and speed, were reported by [Seron et al., 1990; Zienkiewicz and 
Taylor, 1991; Bathe, 1990;]. It was concluded that, even though central difference 
scheme is the best in terms of computational speed and cost, Newmark’s constant 
average scheme gives better accuracy and offers unconditional stability. 


1.3 PRESENT WORK 

The present work deals with study of crack in plates when geometric non- 
linearity is considered under a static load and dynamic load (Suddenly applied load). 
A code developed by Dubey[1993] for linear plate analysis is taken as basis to further 
implement geometric non-linear analysis. The present work makes use of a mixed 
Finite Element Analysis using Reissner-Mindlin Plate theory is used to model the plate 
bending. To avoid locking or singularity a robust triangular element (T6/3 B3) 
[Zienkiewicz and Taylor, 1991] is used. Newmark’s implicit scheme is used for time 
integration. Chapter 2 deals with theoretical part of the formulation of the plate 
bending, geometric non-linearity and the iterative scheme employed. Chapter 3 
describes the Finite Element implementation of geometric non-linear analysis and 
integration with dynamics. Chapter 4 deals with the validation of the code and 
behaviour of plates under dynamic loading when geometric non-linearity is considered. 
Analysis of plate with crack under dynamic loading was also discussed. Chapter 5 
deals with the conclusions of the present work and the scope of further work. 
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CHAPTER 2 


BASIC EQUATIONS OF GEOMETRIC NON-LINEAR 
DYNAMIC PLATE ANALYSIS 

This chapter presents the basic equations governing the thick piate under 
geometric non-linearity and dynamic conditions. Based on earlier chapter, a 
formulation based on Mindlin plate theory can be used for thick and thin plates. 


2.1 GOVERNING EQUATIONS FOR FLEXURE OF THICK PLATE 

The Mindlin plate theory makes the following assumptions, 

1 . The strain and the stress components in Z-direction are negligible (a^ and 
» 0) Fig. (2.1) where Z-direction represents the thickness direction. 

2. The normals to the middle plane of the plate remains straight during 
deformation. 

Applying above assumptions, equilibrium equations for plate bending problem 
can be obtained [Timoshenko, 1959] in terms of midplane rotations (6^, By), lateral 
displacements (w) of the midplane and the shear forces (S^, Sy) at the cross-section. 
Sign convention and the variables are shown in Fig. (2.1). 

Stress resultant moments and shear forces are defined in terms of stress 


variables as follows: 



(2 . 11 a) 
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A/2 
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a„ Z dz 
A/2 y 
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Z dz 


(2.11 b) 




A/2 

A/2 ^ 



dz 


(2.11 c) 


For the dynamic analysis inertia forces due to mass and moment of inertia are 
considered in equilibrium equations Formulating equilibrium equation using 
d’Alembert’s principle, the following three equations are obtained. 


dx 


+ - s^- p-^ = 0 

BY * ^ 12 * 


(2 . 12a) 


dx 


+ 




0 


(2 . 12b) 


35^ 

dx 


+ 


dr 


+ g- phw 


0 


(2 . 12c) 


where 6 ^, '’ 6 ,,, w represent acceleration of respective variables 0 ^, 6,,, w and 

X y X y 

strain at a height of z from the centre plane are given as, 
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2^ 


dx ' 
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(2 . 13a) 


«y = 


By 


Y yz 


6w 

Bx 


j ^ a . Bw 
fy^ - 0y + -gp 


(2 . 13b) 




ay dx_ 


(2 . 13c) 


2.2 THE MIXED FINITE ELEMENT FORMULATION 

The finite element formulation involves expressing the three sets of variables 8, 
S, w, at any point within the element in terms of nodal parameters of the element using 
shape functions. 


e = 



Ng F 


(2 . 14a) 



(2 . 14b) 


W = N„W 


(2 . 14c) 


where 8, S, w are the nodal values of respective variables. 

Now using eqns. (2.13), (2.11) and eqn. (2.14) can be rewritten as follows: 

F - + w = 0 (2.15a) 
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(2 . 15b) 


L 5 = 0 

12 


+ gr - phv?" = 0 


(2 . 15c) 


Where L, v operators and matrices are given by, 

a . 


L = 


dx 

0 

d 




• 3 ■ 

3 


dx 

dy 

/ V = 

_a_ 

_d_ 


. 9y . 

dx 



C, = 


1 

Gh 


1 0 
0 1 


(2 .16) 


D = 


Eh- 


12(l-v2) 


1 V 0 
V 1 0 


(2.17) 


Now by Galerkin method, using eqns. (2.14) and (2.15), and then integrating over the 
area of the element, the following equations are obtained: 


= Te 


(2.18) 


- PS + Cw = 15 


(2.19) 


+ jFV = 


( 2 . 20 ) 
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Where, 


A = fjLN^)^D{LN,i,^,) dQ , B = f^N^^N.dQ 


C = f^Nf i'7N^) dQ ^ ^ fa 


^0 



f.= 



q dQ 


^ /o'^e'’ ^sdQ , F - p hj^Nj N^ dQ 

In the above equations the vector T represents two moment components 
imposed on the boundary by the traction and q is the intensity of the imposed lateral 
force. The shape functions Nq, are chosen to have continuity, but because of 
shear forces are defined only for inside nodes of the element, Ng can be discontinuous 
between the elements. Such a choice allows S to be defined locally for each element 
and thus be eliminated when P is a non-singular matrix. Finally the problem includes 
only e and w as variables and the equilibrium equations can be expressed in terms 
of matrices as 


A+BP-^B^ BP-^C' 



E 

o' 

Pi 


>0 

C-rp-ijgr C^P-^C 

w_ 

+ 

.0 

F 

% » 

y. 
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w+ ±Ldx 
dx 


Fig 2. 



(a) In-plane and bending resultants for a flat plate. 

(b) Increase of middle surface length due to lateral 


displacemei 



In general eqn.(2.21) can be written as, 


m* [u] + m* [l3] = [ ] (2.22) 

Where u consists of 6^^, 6y, w nodal variables and [K] and [M] are the stiffness and 
mass matrices which can be obtained from equation (2.21). 


2.3 MODIFICATION DUE TO GEOMETRIC NON-LINEARITY 

A plate with u,v in-plane displacements in X and Y directions, and w is deflection 
in Z direction is subjected to in-plane and lateral forces is shown in Fig (2.3a). When 
the displacements are not very small the linear stress strain relations are not sufficient 
in predicting displacements or stresses. The in-plane and lateral deformations can not 
be dealt separately and they are coupled The coupling is due to consideration of 
higher derivative terms. 


From Fig (2.3b) it is seen that displacement w produces some additional 

I 

extension in the x and y directions of the middle surface in addition to in-plane 
extensions u, v and the length dx stretches to 


= 



( 2 . 23 ) 


dbc^ = dx 


h + ±1^? - 

1 / dw\‘ 


8 \ dxj 


+ 


( 2 . 24 ) 
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where ‘w’ is lateral deflection of the plate, 


‘dx’ is length of a strip when the lateral deflection is zero, 

‘dx^’ is the elongated length of the strip after certain lateral deflection. 

Including one higher order terms the total elongation in x direction can be written as 





( 2 . 25 ) 


Where, ‘u’ displacement in X direction. 

First term in the above equation (2.33) is due to in-plane loads, second term is 
due to lateral deflection. The second term in the above equation is taken to be 
negligible in case of small deflections. The second term in the above equation couples 
in-plane and lateral displacements. 

In a similar way higher derivative terms are considered for other components. 
The plate strains are defined in terms of mid surface displacements and curvatures. 


Strains and curvatures are written as , 


ey 

yxy 


[ e ] = 


'xy J 


du 


1 / dw\^ 

dx 


2\dx} 

_dv 


1 (dw 

dy 


2 Uy 1 

du du 


ldw\j 8iv\ 

dy dx 


JUy) 
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_ 

dx^ 


0 

d^w 





u 



0 

dxdy . 




( 2 . 26 ) 
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To be concise above equation can be written as 


e 




( 2 . 27 ) 


Where superscript ‘p’ represents in-plane, ‘b’ represents bending and subscript ’nl’ 
represents non-linear components 

The in-plane stresses and moments resultants are represented in vector form 


as 


[o] 


My 



( 2 . 28 ) 


The ‘stresses’ are defined in terms of the usual resultants as 

Ty = ^yt 

( 2 . 29 ) 


Where ffS:, ay, are the average membrane stresses and, ’t’ is the thickness of the 

xy 

plate 
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Assuming material to behave linearly the stress resultants and strains can be 


related by stiffness matrix D, defined as 


D 


DP 0 
0 p* 


( 2 . 30 ) 


Where the in-plane component of D is 


DP 


[i V 


E 

1-v^ 


V 1 
0 0 


0 

0 

1-v 

2 


( 2 . 31 ) 


and flexure component of D is 




1 V 


Et^ 

12 (l-v^) 


V 

0 


1 

0 


0 

0 

1-v 

2 


( 2 . 32 ) 


The displacements are related to nodal parameters by shape functions, as 


= JWa 


( 2 . 33 ) 
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Where N is shape function matrix which is composed of 


N, 


Nf 0 
0 N^’ 


( 2 . 34 ) 


Where V corresponds to particular node, 
N P is in-piane shape function, 

Nj*^ is bending shape function. 


a 


i ~ 



( 2 . 35 ) 


Where a P is nodal parameters of in-plane, 
a,^ IS nodal parameters of bending. 
Which can be represented as 


. ^i. 

g 

II 

Ki' 






( 2 . 36 ) 


e^, 0„ are the rotations in x and y directions. 

X y 

Wj is lateral deflection. 
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2.4 FINITE ELEMENT PROCEDURE FOR INCLUSION OF GEOMETRIC NON- 
LINEARITY 

The finite element procedure for non-linear analysis is based on the fact that the 
displacements are large and equilibrium conditions between internal and external 
‘forces’ have to be satisfied 

Hr(a) = f a dV - 1 (2.37) 

J V 


t|r sum of external and internal generalized forces 
B is defined from strain displacement definition as 

de = jB da (2.38) 

B = + B^J (a) (2.39) 

Bq is due to small deflection theory (that is due to first term of equation (2.33), 
Bf^i is due to geometric non-linearity effect (that is due to second term in equation 
(2.33) ). 

Thus by taking appropriate variations with respect to a of equation (2.41) it can 
be written as 

=( dB o dv + do dv = Kj. da (2.40) 
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Stress is related to strain by relation 


a = D e , do - -Dde 

By using above relations di|r is written as 

= j dB^x o dv + K da 

Where 

K = fT^DBdv = + JCjjj 

KqP in-plane stiffness matrix 
bending stiffness matrix 

= f (bJ d b^x * d * bZ. d dv) 

J V 


f dB'^ o dv = K. da 
J V 

I 


cfi}r = (JC^ + + K^i)da = da 

Kpi stiffness matrix due to geometric non-linearity. 

Kj is total are tangential stiffness matrix 

Ka is stiffness matrix dependent on stress level. 


centum. 


( 2 - 41 ) 


( 2 . 42 ) 


( 2 . 44 ) 


( 2 . 44 ) 


( 2 . 45 ) 


( 2 . 46 ) 





The above procedure for calculating tangential stiffness matrix in geometric non 
linear problem can be put in the following algorithm. 

(1) Evaluation of the linear small deformation matrix. 




Ki 0 


0 kL 


(2.47) 


Bi dv 

J V 


:2 .48) 


KP is in-plane stiffness matrix. 

is bending stiffness matrix calculated from equation (2. 28) 

(2) Evaluation of higher order deformation stiffness matrix 

= f (^o"' B n B Bj dv (2.49) 

J V 


where B^ 


Bi 

0 

il 

0 

bA 

0 

B^\ 


.0 

0 . 


(2.50) 


On substituting and in above equation can be written as 





0 


sym. , 




dv 


(2.51) 


B^i is due to higher order deformation . 
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(b) Evaluation of B^, 

From the strain-displacement relations 


du ^ 


zz. ^ — 

^ non- linear ' 2 ^ j 

dy 2lay] 

_ _1 / dwW 

$y 9x 2 \3xA3y/ 


non-linear term can be v\/ritten as 
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Where for the sake of simplicity it is assumed that 
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0.^ = 
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L®y J 






= AG 


(2.52) 


(2.53) 


(2.54) 


(2.55) 


(2.56) 
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(c) Evaluation of B^p 


B 


P 

O 


d 

dx 

0 


0 


A 

By 


0 , 


_0_ 

By Bx 


0 , 


N^, 0, 

0/ 


(2.57) 


( 3 ) 


Where {N.,,N 2 are shape functions of and By respectively. 

Evaluation of 

is dependent on the stress level, which found from the following equation 


iC 


- L 


0 0 
dA^ 0 


■^xy 


M, 


^ J 


dv 


(2.58) 


After some manipulation of the above equation it can be written as 






T T 

•‘‘X ’*‘xy 

T T 


G dv 


(2.59) 


Where G is a shape function matrix of 6^, and Gy written as 

'n, 0 0 0 . . / 

(? = ^ ^ ^ 

0 0 Wj 0 W3 , . . 


(4) Evaluation of tangential stiffness matrix 

After evaluating in-plane and bending and stress dependent stiffness matrices, 
the total tangential stiffness matrix is calculated by this equation 

JCj. = + (2.61) 
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This tangential stiffness matrix is used in finding non-linear displacements in 
plates. The next session deals with the iterative scheme used for getting non-linear 
displacements. 

2.5 ITERATIVE SCHEME 


Iterative scheme employs modified Newton Raphson method. In this method 
variable stiffness matrix is replaced by a constant stiffness matrix for a given load step. 
Even though convergence is slower in this method it reduces the computational time 
needed. 

The basic equation for doing this iterative method 

[U]i = [U (2.62) 

Left hand side of the eqn. (2.62) represents displacement after the i^^ load step, 
and right hand side is sum of total displacement after (i-1)^^ load step and 
displacement due to i^*^ load step. 

Kg = f{U) (2.63) 


Where U represents displacement 
Kq represents global stiffness matrix. 

du^= m2\[dF-dfJ (2.64) 

dF is the incremental applied load. 

df^ is the resultant load due to internal stresses after m equilibrium iterations, 
m is equilibrium iteration number. 
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The stiffness matrix used in the above equation is stiffness matrix corresponding 
to (i-1)*^ load step. But for calculating unbalanced loads it is required to calculate the 
stiffness matrix in each equilibrium iteration, which consumes considerable time. In 
order to reduce computational time unbalanced loads are calculated at elemental level 
and assembled. 

The above eqn. (2.64) can be written as 

In the above eqn. (2.65), residual forces are calculated for every element and 
assembled. In the above equation [Kq]^.., is updated stiffness matrix after (m-1) 
equilibrium iterations. 

convergence criterion : 

The error in each iteration is represented as follows, 

e = [dcg - [dU^J 


Convergence is said to be achieved when 



^ tt 


Where a is predefined convergence limit. 
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CHAPTER 3 


IMPLEMENTATION OF FINITE ELEMENT CODE 

This chapter explains the different aspects of numerical method in developing 
the FEM code, including element selection. 

3.1 ELEMENT SELECTION 

It is well known that complex geometric domains can be easily discretized into 
triangular elements. In literature for different types of problems, different types of 
elements are suggested as they should be free from defects of locking and 
singularity.The performance of the element under different type of loading and 
boundary conditions should also be stabIe."Robust triangular element" with suitable 
nodal variables satisfy the stability conditions. [DUBEY,1993] 

The choice of the shape functions and variables should be such that the 
Babushka-Brezzi[Zienkiewicz and Lefebvre,1988] conditions are satisfied by the 
system and thus, stability is ensured.The most vital and strictly necessary condition for 
the non-singularity of the system given by 




P 



( 3 . 2 ) 


Where ng, ng, and n^,, standfor the number of variables for rotation, shear force and 
lateral displacements respectively in an element. If either of the two conditions are 
violated, the locking behaviour will occur in the solution. 
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For the dynamic analysis if an implicit time integration scheme is used, it is 
necessary to use higher order finite elements and a consistent mass matrix. The 
higher order elements are effective in the representation of bending behaviour.These 
elements should be employed with consistent load vector, in that the mid side and 
corner nodes are subjected to their appropriate load contribution in the analysis. 

In the mixed finite element formulation, triangular element (T6/3) has six 
boundary nodes which are used to determine a quadratic variation of 6 and w 
[Zienkiewicz and Lefebvre,1988]. Shear forces are determined by a linear field by 
three internal nodes with no element inter connections fig(3.1a). 

When standard patch tests [Zienkiewicz and Lefebvre,1988;Zienkiewicz et 
al.,1986] are used with this element, the results indicate that, the element has incipient 
locking possibility andthus it is not a satisfactory element. But the examination of 
results indicate that the satisfaction of the requirements can be achieved by the 
addition of the 8 variables in the interior of the element. The three internal nodes 
created have the shape functions used of the form L.,^L 2 L 3 , provide a new element 
(T6/3 B3) which passed the same patch test successfully in all cases [Zienkiewicz and 
Lefebvre,1988], shown in fig(3.1b). 

I 

This particular element is chosen for the present analysis. It has 30 degree of 
freedom, and after the condensation of shear force terms the number of variables 
reduces to 24. And for consideration of in-plane displacement two degree of freedom 
are added to each of the 9 nodes, thus increasing total degree of freedom to 42, 
shown in the fig (3. 1c) 
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3.2 TIME INTEGRATION SCHEME 


The equation (2. 29) is integrated by a numerical procedure. The integration is 
based on two ideas. 

1 . Static equilibrium including inertia force is achieved at discrete time intervals At. 

[jq t+At|-y] ^ [j;^t+At (3.3) 

2. Variation of displacement, velocity and acceleration within each time step is 
assumed. This assumption decides the accuracy and stability of the solution 
procedure. 

The idea no 2 decides the time integration scheme. The choice of a scheme 
depends on the finite element idealisation. However, the finite element idealisation to 
be chosen depends in turn on the actual wave velocity of the medium to be analyzed. 
It follows, therefore the selection of an appropriate finite element idealization of a 
problem and the choice of effective integration scheme for the response solution are 
closely related and must be considered together. 

Out of many schemes ’Newmark’s implicit scheme is most accurate and 
unconditionally stable. 

3.3 NEWMARK IMPLICIT SCHEME 

The Newmark’s constant acceleration scheme is written as follows [Bathe, 1 990] , 

^ [ft + [(1-6) + 6 ] At (3.4a) 


^ |- (A-oc) 1^^= + a ] At^ (3.5b) 

2 
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Where a and 6 are parameters that control integration accuracy and stability. 
Newmark originally proposed an unconditional stable scheme, in which case 6 = 0.5 
and a = 0 25. Using eqn(3.3) and eqn(3.5) final form can be obtained as, 

[if + a, iy] + im {a,u^ ^ a, u^) (3 .ea) 


- u^) - - a^u^ (3.6b) 


u 


t + At 


+ a. + a^ u 


t-i-At 


(3.6c) 


^0 




(3.7) 


The implicit scheme, as can be seen from the final formulation requires 
inversion of [K+ag M] matrix. For this purpose Cholesky decomposition technique is 
computationally efficient. 


3.31 Time step 

After deciding on time integration scheme, it is necessary to use a time step 
suitable with finite element mesh. If the critical wavelength in the medium is denoted 
by L^, then the time step can be decided as , 
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And the effective length of a finite element should be 


L, = C At 


Where ’c’ is wave speed and ’n’ is the number of nodes per critical wave length. 
The effective wave length and corresponding time step must be able to represent the 
complete bending wave accurately and is chosen appropriately depending on the kind 
of element idealization and time integration scheme used. It is reported that at least 
8 nodes per shortest wave length are required in order to produce all artifacts of wave 
propagation [Seron, 1990] 

However, for higher order (parabolic and cubic) continuum elements the time 
step has to be further reduced, because the interior nodes are stiffer than the corner 
nodes [bathe, 1990]. 

3.4 COMPUTER IMPLEMENTATION 

The computer code developed by Dindore[1994] has been enhanced to take 
into consideration of geometric non-linearity in the present work. Since geometric non- 
linearity is considered, load applied on plate should be in small increments and for 
each increment of load, convergence criteria has to be checked. In dynamic loading 
the time step used should be smaller when geometric non-linearity is considered than 
without, since the structure stiffens with the displacements. 

In case of crack in plate under impact loading, it is assumed that an initial crack is 
present and it is extended and restored at every time step. Energy before crack 
extension and after crack extension is calculated, the difference of these two energies 
is the energy release. 
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Energy in plate is calculated by the expression, 

p.E= ^[uv miu] ^ - lpieu ]^ o . s ) 

Where [u] is velocity vector, P is load vector, [K] is the stiffness matrix, [M] is the 
mass matrix. 
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CHAPTER 4 


RESULTS AND DISCUSSIONS 


This chapter is divided into two parts, first consists of validation of the present 
method while second part analyses the effect of geometric non-linearity in plate, with 
and without crack. 

4.1 VALIDATION OF THE METHOD 

Validation of the present method is done in following cases. 

Case 1 • Stresses and deflections for a elastic problem are calculated for standard 
cases (Timeshenko, 1959). A problem of fixed plate subjected to concentrated load 
is considered (fig 4.1). 

The material and geometric properties are as follows: 


Modulus of elasticity. 

E 

= 200 Gpa 

Poisson’s ratio. 

u 

= 0.25 

Density, 

P 

= 7800 Kg/m® 

Thickness of the plate. 

h 

= 0.02 m 

Length of the plate, 

L 

= 2 m 


Considering the symmetry of the problem, only one quarter of the plate is analyzed. 
Mesh used is shown in the Fig 4.2. The mesh used is generated and optimized by 
an interactive preprocessor already developed (Kishore, 1996). An enlarged view^ 
showing details of the quarter plate is shown in Fig 4.3. Compared to the theoretical! 

I 

value of central deflection of 1 .52cm for a load of 0.1 MN, the F.E.M result obtained is^ 
1 .50cm. As can be seen, the value obtained by F.E.M is within 1 .5% of the theoretical 
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value and is reasonably accurate. 


Case 2: Deflections are calculated by considering geometric non-linearity effect for a 
standard case of rectangular plate fixed on all sides and subjected to uniform load. 
In this problem mesh used is same as that of in fig(4.2) and the material properties 
used as in case 1 . Fig 4.4a shows the theoretical results given in Timoshenko and 
Krieger[1989]. Fig 4.4b shows the results obtained by the present method. Compared 
with the theoretical value of the central deflection, /h, of 1.52 the F.E.M results 
obtained a value of 1.49 for the same load value, qbVDh, of 150. As can be seen the 
error of the present method is below 2%. 

Case 3: The response of an elastic cantilever under dynamic condition was determined 
by the present method and compared with those given by Bathe and Ramm [1 975] 
for the purpose of validation. 

The material and geometric properties used are as follows: 

Modulus of elasticity, E = 8.277 Mpa 

Poisson’s ratio, o =0.2 

Density, p =10.7 Kg/m^ 

Thickness of the plate, h = 0.025 m 

Length of the plate, L = 0.25 m 

Uniformly distributed load, Qq = 20292 N/m^ 

The load is a suddenly applied step load shown in Fig 4.5a. 
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The variation of dynamic load is assumed as follows. 


For load #1 , 


Q = 

0 

1 

o 

at 

t 

= 0 

<7 = 

<7o 

at 

t 

fV 

O 


(4.1) 


This loading as shown in fig(4.5a) 


For load #2 


<7 = 


1 


t^RT 

<7 = 

2 

<7o 


\RTl\ 

t > RT 


(4.2) 


where RT is the rising time. 

Mesh used is same as in Fig 4.3. 

The element size and the time steps are chosen in such a way so that the 
solution remains stable. Time step was selected as a fraction of fundamental time 
period. In this, time step is taken as 45jiS. As explained in earlier chapter the results 
were obtained using Newmark integration scheme. 

Fig (4.6) presents the end deflection of the cantilever. For the present method 
gives maximum deflection, of 0.76 which compares well with the value of 

0.745 obtained by Bathe and Ramm [1975]. 

4.3 EFFECT OF GEOMETRIC NON-LINEARITY IN PLATE UNDER IMPACT 
LOADING 

Fig 4.7 presents the end deflection of a cantilever subjected to impact loading 
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load #1 shown in Fig. 4.5a with geometric non-linearity. The figure also shows the 
results obtained for a plate without geometric non-linearity. As expected, it is seen in 
this figure that the cantilever stiffens with increasing displacements and there is 
decrease in amplitude and effective period of vibration. 

When higher impact are loads applied the effect of geometric non-linearity is 
increased as seen in figures (4.7b) and (4.7c). 

4.4 GEOMETRIC NON-LINEARITY EFFECT ON CRACK IN PLATE UNDER 
IMPACT LOADING 

Fixed plate with a crack subjected to uniformly distributed impact load is shown 
in Fig 4.12. The plate is analyzed for two load cases: 

Case 1 : In this impact load #2 as shown in Fig 4.5b was applied. In this problem it was 
assumed that an initial central crack of length 71.4cm of the edge length was 
assumed. Time step of lOjis was chosen. Mesh used was same as in Fig 4.2. 

The material and geometric properties are as follows: 

Modulus of elasticity, E = 200 Gpa 
Poisson’s ratio, u = 0.25 

Density, p = 7800 Kg/m^ 

Thickness of the plate, h = 0.02 m 

Length of the plate, L = 2 m 

Fig 4.9 gives the variation of energy release rate as a function of time, t v® 
observed that energy release rate is higher with geometric non-linear consideration 
than without geometric non linear consideration by 10.7 percent. This is due to the 
contribution of energy due to membrane stresses. These membrane stresses are 
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neglected in linear problem. 


Case 2: For the same material properties, mesh and crack conditions as above and 
an impact load given by load #1 (Fig 4.5) the plate was analyzed. Fig 4.10 shows the 
variation of energy release rate as a function of time. It is observed in that energy 
release rate in both linear and non-linear case is high due to than case(1), in linear 
case it is by 5% in non-linear case it was by 10%. This is due to the extra energy that 
has been put in due to sudden initial rise of the applied load. The energy release to 
is higher with geometric non-linear consideration than without geometric non-linear 
consideration by nearly 18 percent. 
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Plate subjected to a concentrated load p at the 
center . 


Figure 4.1. 




Figure 4.2. Mesh of 81 elements on the quarter portion of plate. 




Figure 4.3. Enlarged quarter portion of the plate. 
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Figure 4.7a Finite displacement dynamic response of a 
cantilever under uniformly distributed load 
q = 0 . Iq^ , with geometric non-linearity 
consi derati on . 
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Figure 4.7c Finite displacement dynamic response of a 

cantilever under uniformly distributed load 
q= 0.?S with aeometric non-linearity 
consideration. 
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Figure 4.9 Variation of energy release rate as a function of 
time for a clamped plate under dynamic loading(Load t2). 
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Figure 4.10 Variation of energy release rate as a function of 
time for a clamped plate under dynamic loading(Load #1). 



CHAPTER 5 


CONCLUSIONS AND FUTURE WORK 

Based on the present finite element analysis the effect of geometric non-linearity 
on the behaviour of a plate with and without crack subjected to sudden loading followi 
ng conclusions can be drawn; 

1 . Static and dynamic deflections of plate with geometric non-linearity are lower 
than the plate without geometric non-linearity consideration. 

2. The frequency of oscillations of plate with geometric non-linearity are more than 
plate without geometric non-linearity consideration. 

3. The effect of geometric non-linearity is highly dependent on the amount of 
displacements. 

4. Energy release rate is 10-20 percent higher when geometric non-linearity is 
considered than, linear case under impact load conditions. Thus, it is important to 
include this effect while investigating crack initiaion. 


SCOPE OF FUTURE WORK 

The mesh can be made finer and the time step can be made smaller to observe 
the sharper impact loading on plates with geometric non-linearity. It can be further 
modified to include large deformation and can be integrated with material non-linearity. 
This analysis can be tried to various configurations such as spherical shells and 
arches. 
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